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ABSTRACT 

Given a set of astrometric observations assumed to belong to the same ob- 
ject, the problem of orbit determination is to compute the orbit with all the 
necessary tools to assess its uncertainty and reliability. Under the conditions 
of the next generation surveys, with much larger number density of observed 
objects, new algorithms, or at least substantial revisions of the classical ones, 
are needed. The problem has three main steps, preliminary orbit, least squares 
orbit, and quality control. The classical theory of preliminary orbit algorithms 
was incomplete, in that the consequences of the topocentric correction had not 
been fully studied. We show that it is possible to rigorously account for the 
topocentric correction, possibly with an increase in the number of alternate 
preliminary orbit solutions, without impairing the overall orbit determination 
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performance. We have developed modified least squares orbit determination 
algorithms, including fitting methods with a reduced number of parameters 
(required when the observed arcs have small curvature), that can be used to 
improve the reliability of the orbit computing procedure. This requires suitable 
control logic to pipeline the different algorithms which we have defined and val- 
idated through numerical simulations. We have tested the complete procedure 
on two simulations with number densities comparable to that expected from 
the next generation all-sky surveys such as Pan-STARRS and LSST. To con- 
trol the problem of false identification (where observations of different objects 
are incorrectly linked together) we have introduced a quality control on the fit 
residuals based upon an array of metrics and a procedure of normalization to 
remove duplications and contradictions in the output. The results confirm that 
large sets of discoveries can be obtained with good quality orbits and very high 
success rate losing only 0.6 to 1.3% of objects and a false identification rate in 
the range 0.02 to 0.06%. 

Key Words: Celestial Mechanics; Asteroids, Dynamics; Orbits 

1 The Problem 

The problem of preliminary orbit detcrminatior{3 is old, with very effective so- 
lutions developed by [Laplace 1780| and [Gauss 1809j . Of course the methods 
of observing Solar System bodies have changed radically since classical times 
and have been changing even faster recently due to advances in digital astrom- 
etry. The question is, what needs to be improved in the classical algorithms to 
handle the expected rate of data from the next generation of all-sky surveys? 
Alternatively, what can we now use in place of the classical algorithms? 

The issue is not one of computational resources because these grow at the 
same rate as the capability of generating astrometric datgd. Reliability is the 
main problem when handling large astronomical data sets (millions of individual 
detections of Solar System objects). An algorithm failing once when used 1, 000 
times may have been considered perfectly reliable only a few years ago but in 
the present situation we must demand better performance, and even more so in 
the near future. 

This is particularly important because of the strong correlation between 
difficulties in the orbit computation and the scientific value of the discovered 
object. Main Belt Asteroids (MBA) are commonplace and their orbits are easily 
computed. Only a few in a 1, 000 of the objects (to a given limiting magnitude) 
are the more interesting Near Earth Objects (NEO) while few in 100 are the 
equally interesting Trans Neptunian Objects (TNO); in both cases the computa- 
tion may be much more difficult for reasons explained later. Thus an algorithm 

'^Also called Initial Orbit Determination (lOD). 

^Moore's empirical law predicts an exponential growth of the number of elements on a 
chip with time and this affects the number of pixels in a CCD and the performance of the 
computers used to process astrometric data in the same way. 
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computing orbits for 99% of the discoveries may be failing on a large fraction of 
the more interesting objects like the NEO and TNOs. 

To find reliable algorithms we need a good qualitative understanding of the 
solutions and a firm control of the approximations used. The classical algorithms 
contain several approximations, we will show that one of them is the most 
dangerous: neglecting the topocentric correction and assuming that the observer 
sits at the center of mass of the Earth. The well known method of reducing to the 
geocenter by assuming the object's distance and then iterating the preliminary 
orbit computation does not solve the problem. 

Although Laplace's and Gauss' methods each have their supporters it is gen- 
erally believed that they are equivalent to a good approximation. In fact, we 
show that this is true only when the topocentric correction is neglected. As 
already pointed out by [Marsden 1985] . Gauss' method has significant advan- 
tages in the way it can include the topocentric correction. However, the classical 
qualitative theory |Charlier 1910] of the number of alternate solutions for the 
preliminary orbit applies only to Laplace's method neglecting the topocentric 
correction. 

We show that a reliable method, effective under the present observing con- 
ditions, can account for the topocentric correction in Gauss' method although 
such a correction could also be included in Laplace's method. Then we need a 
qualitative theory, replacing the one of Charlier, for Gauss' method that does 
not assume geocentric observations. We have developed such a theory and found 
that the number of alternate preliminary orbit solutions can be larger than in 
Charlier's theory, namely, there may be double solutions at opposition and triple 
solutions at low solar elongations. 

A reliable orbit determination algorithm should begin with a preliminary 
orbit algorithm that accounts for the possibility of double and even triple solu- 
tions. This is especially important to reliably handle NEO discoveries which are 
affected by non-unique solutions (occurring in most cases near quadrature and 
sometimes even at opposition) . Moreover, particular provisions are required for 
the case in which the quantities used as input to the preliminary orbits, the 
curvature components, are poorly determined to the point that their sign is un- 
certain. This happens when the observed arc is too short and when the observed 
object is very distant. For the case in which TNOs are being discovered, weakly 
determined or totally undetermined preliminary orbits are the rule. Two classes 
of methods can be used in such cases: the Virtual Asteroids (VA) methods and 
the constrained least squares solutions. 

A Virtual Asteroid is a fully specified orbit with six orbital elements that 
is compatible with the available observations but by no means determined by 
them. In different types of VA methods one to thousands (depending upon the 
purpose) of VA are selected either at random or by some geometric construction 
inside the region in the space of the orbits compatible with the observation^. 
There are more than half a dozen VA methods available in the literature and we 

3 The method of IVaisala 19411 is one way of selecting just one orbit when there are only 
two observations. The modern VA methods are generalizations of the one by Vaisala with the 
advantage that they also perform well far from opposition. 
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will not review all of thenfl but just present one version which is particularly 
effective for the problem we are discussing in this paper. It consists in choosing 
just one VA in such a way that both the condition of being compatible with 
the observations and having an elliptic orbit are satisfied. It is derived from the 
theory of the Admissible Region we have developed [Milani et al. 2004] . This 
is especially effective when the preliminary orbits computed with the classical 
algorithms (or even modern versions) arc ineffective as happens in most cases 
for TNOs observed near quadrature. 

The main purpose of preliminary orbits is to be used as a first guess for the 
nonlinear optimization procedure [Differential Corrections) that identifies the 
nominal orbit fulfilling the principle of Least Squares (of the residuals). If the 
preliminary orbits are well defined by the curvature of the observed path, then 
one of them is likely to be close enough to some least squares orbit to belong to 
the convergence domain of the differential corrections; then it is sufficient to use 
all the preliminary orbits as first guesses. If the preliminary orbits are poorly 
defined they are likely to be so "wrong" that they will not lead to convergence 
of the differential corrections, an inherently unstable procedure when using a 
starting point far from the nomina{f|. 

Thus it is essential to increase the size of the convergence domain by using 
modified differential corrections methods. Many of these methods exist and 
most of them have one feature in common: the number of parameters determined 
is less than 6. There are 4- fit methods in which 2 variables are kept fixed (at 
the value determined by the preliminary orbit or by some other criterion) and 
the others are corrected in an iteration converging to the minimum of the sum 
of squares of the residuals restricted to a 4-dimensional submanifold of the 6- 
dimensional space. There are 5-fit methods in which one parameter is fixed 
although it does not need to be one of the orbital elements but might be defined 
in some intrinsic way adapted to the particular problem at hand. The 5-fit 
method we use is fully documented in [Milani et al. 200"5a) . 

Although the focus of this paper is on new theory and on algorithm doc- 
umentation we felt the need to thoroughly test the actual performance when 
our methods and software are confronted with the expected data rate of the 
next generation surveys. Thanks to the collaboration with the Pan-STARRS 
[Jedicke et al. 2007] and LSST [Ivezic et al. 2007j projects we have had the op- 
portunity to use simulations of future survey observations to precisely measure 
performances using the simulation source catalog as ground truth. In particular 
we report here on two tests, one based on a small but focused simulation con- 
taining only the most difficult orbits (NEOs and TNOs), and one representing 
the full data rate from a next generation survey containing all classes of solar 
system objects (of course the numbers are dominated by MBA). We believe the 
results are a convincing confirmation that we are ready to handle the data from 
the next generation surveys. 

*For a recent review see jMilani 2005) . 

^This property is shared by all variants of Newton's iterative method that is notoriously 
prone to chaotic behavior when used outside the convergence domain. 
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2 Equations from the Classical Theory 



There are so many different versions of preliminary orbit determination methods 
and there is so httle in the way of a standard notation that it is not possible 
to simply copy the equations from some reference. Thus, we have chosen to 
provide here a compact summary of the basic formulae, in particular discussing 
the dynamical equation and the associated polynomial equation of degree 8 for 
both methods developed by Gauss and Laplace. We also summarize Charlier's 
theory on the number of solutions for Laplace's method, to be compared to the 
new qualitative theory, including topocentric correction, of Section [J] 



2.1 Laplace's Method 

The observation defines the unit vector p — (cos (5 cos a, cos (5 sin a, sin (5) where 
(a, S) are the topocentric right ascension and declination. The heliocentric po- 
sition of the observed body is 

r = p + q = pp + qq 

where q is the observer's positior0. Let s be the arc length parameter for the 
path described by the relative position p{t) and 77 the proper motion 



ds r~ ~ ^ d Id 

-77 ya cos^ 6 + 6^ ; — ^ -— . 

at ^ ds rj dt 

We use the moving orthonormal frame [Danby 1962[ Sec. 7.1] 

. , dp ^ ^ , 
p,v=— ,n = pxv (1) 

and define the geodesic curvature n by the equation 

dy' . . 

- = -p + .n. (2) 

Then the relative acceleration is 
d'^p 

^ = (P - pjf)p + (pi + 2p?7)v + {pifn)h (3) 

and the differential equations of relative motions are 

d^p .. .. UQ ur 

dt'' q-^ r"^ 



(4) 



with the following approximations: q = q® coincides with the center of mass of 
the Earth, the only force operating on both the Earth and the object at r is the 



"Aberration, a consequence of the finite speed of light c, is accounted for by assuming a 
time at the object t^i^j = t ~ p/c different from the observation time t. 
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gravitational attraction by the Sun, i.e. there are no planetary perturbations 
(not even indirect perturbation by the Earth itself). From (jSj-n =([4])-n 

d^P „2 , .. .... n 1 



which can be presented in the form 

C^ = l-4 with C=^^ (5) 

referred to as the dynamical equation in the literature on preliminary orbitfl 
C is a non-dimensional quantity that can be when r = q or undetermined (of 
the form 0/0) in the case that the O(At^) approximation fails, i.e. the object 
is on an inflection point with tangent pointing to the Sun. 

Given p, to complete the initial conditions p is solved from ([3])-v =(j4])-v 

- p—^ + p-^ ^ pi] + 2pr] . (6) 



2.2 Charlier's Theory 

Eq. ([5]) is the basic formula for Laplace's method using the solution in terms of 
either p or r which are not independent quantities. From the triangle formed 
by the vectors q, p, r we have the geometric equation 

= p^ + 2pqcose + q^ (7) 

where cose — q • p is fixed by the observation direction (e ~ 180°— solar 
elongation). The level curve C = is the zero circle r = q. By substituting p 
solved from eq. ([71) in ([5|), removing the square root by squaring and multiplying 
by (C 7^ 0, otherwise r q) we obtain the polynomial equation 

P{r) = C^r^ - 9^6(1 + 2Ccose + C^) + 2q'^r^{l + Ccose) - 9^ = . (8) 

Since P{q) = there is the trivial root r ~ q, due to a singularity in the 
spherical coordinates. There can be other spurious solutions of the polynomial 
equation ([8]) corresponding to p < in eq. ([5]) . 

The qualitative theory of [Charlier 1910| on the number of solutions is ob- 
tained by analyzing these equations with elementary methods. The sign of the 
coefficients of eq. ^ is known: —(1 -I- 2Ccos£ -I- C^) < and (1 -I- Ccose) > 
(see [Plummer 1918j ). Thus there are 3 changes of sign in the sequence of 
coefficients and < 3 positive real roots. By extracting the factor (r — q) 

P{r) = {r - q) Pi{r) ; Pi{0) = q' ; Pi (q) = C (C - 3 cose) . 

The number of solutions of the polynomial equation changes where Pi (q) changes 
sign, at C = ^ r = q and at C — 3 cose ~ 0. The latter condition 

^It is, in fact, the component of the dynamical equations along the normal to the path. 
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defines the limiting curve; in heliocentric polar coordinates {r,(f>), by using 
^ + — 2r q cos <j> 

3 

4- 3-COS0 = % . (9) 
q r'^ 




Figure 1: Level curves of C(r, p) (solid lines), limiting curve (labeled), zero circle 
(dashed). For a given value of C and an observation direction (dotted) there 
can be either 1 or 2 solutions, e.g. for C = 0.3 there are 2. 

Following [Charlier 19101 ICharher 19U\ and [Plummer 19T8] the number of 
solutions can be understood with the help of a plot of the level curves of C(r, p), 
in a plane with the Sun at (0,0), the Earth at (g,0) and the position in each 
half-plane defined by the bipolar coordinates (r, p) . The limiting curve and 
the zero circle can be used to deduce the number of solutions occurring at the 
discovery of an object located at any point of the plancH. There is only one 
solution on the right of the unlimited branches of the limiting curve, around 
opposition. There are two solutions for every point of the region between the 
unlimited branches and the zero circle. Inside the zero circle and outside the 
loop of the limiting curve there is only one solution. Inside that loop there are 
always two solutions. 

Note that the classical theory by Charlier assumes that there is always at 
least one preliminary orbit solution. This results from two implicit assumptions: 
that the observed object exists (not being the result of a false identification) and 

*This plane does not correspond to a physical plane in that it also describes the points 
outside the ecliptic plane. 
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that the value of C is measured exactly, or at least to good accuracy, from the 
observations. C contains k which can be difficult to measure from a short 
observed arc, thus both assumptions may fail as discussed in Section [51 and 



2.3 Gauss' Method 

The method by Gauss uses 3 observations corresponding to heliocentric positions 

r. = p, + q. 1 = 1,2,3 (10) 

at times ti < t2 < h with ti — tj = 0{At) ^ period and the condition of 
coplanarity: 

Airi - ra + Asrg = . (11) 

From pTjl xr, • c, where c = x r^, the coefficients Ai, A3 are obtained as 
triangle ratios 



A, = 



r2 X rs • c 



A, = 



ri X r3 • c 
From Uni) and pi x pv pTj) : 



ri X r2 • c 
ri X ra • c 



P2[Pl X P3 • P2] = Pi X P3 • [AiQi - q2 + A3q3]. 



(12) 



Next, the differences — r2 are expanded in powers of tij = ti — tj = 0{At). 
e.g. by using the f,g series formalism = /ir2 + gir2 and Taylor expansions 



/, = l-§| + 0(At3) 



9i 



t,7 1 



6 rl 



-0{At^) 



Then r,; x r2 = -giC, ri x r3 = {figs - hgi)c and 

93 ^ n \ 



Ai = 



/i53 - /351 



> ; A3 = 



/l<73 -/35l=t31 (1- 



/i33 - hai 
0{At^) 



> 



Using ^ and USD in HID 

Ai 

A3 



^32 
t^ 


n 


^^ 


t21 

hi 


n 





©(At^) 
©(At^) 



(13) 

(14) 
(15) 

(16) 
(17) 



Let P be 3 X volume of the pyramid with vertices q, r 1 , r2 , ra P = pi x p2- P3\ 
by substituting it and (fTH]) . (|17p in ffl2\i . with simple manipulations of the timecl 



Pp2t31 = Pi X P3 • (i32qi - <3iq2 + i2iq3) 



(18) 



^Use t|j - = i2i (iai + ts2 ) and - = <32 (tsi + t2i ) 
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^^[^32^21(^31 +t32)qi +^32^21(^31 +i2i)q3] 



If the terms 0{At'^) arc neglected and wc let _B(qi,q3) represent the coefficient 
of the 1 /rf term in (jlSp then 



-B(qi,q3) = 77^32*21/01 X P3 • [(^31 + t32)qi + (*31 + i2l)q3] 

6 



(19) 



Then multiply p8)) by q2/-^(qi:q3) to obtain 



P P2 tsi 3 ^ g| 

S(qi,q3)'^' 



^(qi,q2,q3) 



where 
Let 

and then 



-S(qi,q3) 

^(qi,q2,q3) ^ql pi^ ps- [teqi -i3iq2 +t2iq3]- 

Phiqt J, -4(qi,q2,q3) 



(20) 



Co 



-B(qi,q3) 



^(qi,q3) 



Oo — — "0 T 



(21) 

92 n 

is the dynamical equation of Gauss' method, similar (but not identical) to eq. (O 
of Laplace's method. Using ([7|) at time t2 (with 52, P2, '"2 and £2)'- 

Po{r) = Cy2-ql4ihl + 2CohoCose2 + C^) + 2qlrl{ho+Cocose2)-ql = (22) 

where the sign of the coefficients is as for ([8]) , apart from ho + Cq cos £2 whose 
sign may change depending upon Hq. Note that Po{q) 0, no root can be found 
analytically. The number of positive roots is still < 3 but a qualitative theory 
such as the one of Section 12.21 is not available in the literature. 

After the possible values for r2 have been found the corresponding p2 values 
are obtained from eq. (|2ip and the velocity r2 can be computed, 6.17. from the 
classical formulae by Gibbs [Herrick 19711 Chap. 8]. 



3 Topocentric Gauss-Laplace Methods 

The critical difference between the methods of Gauss and Laplace is the fol- 
lowing. Gauss uses a truncation (to order 0{At'^)) in the motion r{t) of the 
asteroid but the positions of the observer (be it coincident with the center of 
the Earth or not) are used in their exact values. Laplace uses a truncation to 
the same order of the relative motion p{t), thus implicitly approximating the 
motion of the observer. In this section we examine the consequences of the 
difference between the techniques. 



9 



3.1 Gauss-Laplace equivalence 

To directly compare the two methods let us introduce in Gauss' method the 
same approximation to order 0{AP) in the motion of the Earth which is still 
assumed to coincide with the observer. The /, g series for Earth are 



fJ-th [ 3(q2 ■ q2)q2 



9| 



q2 



(23) 



By using ^ in (dH) we find that 



^(qi,q3) = ^i32t2lPl X P3 ■ [3t3iq2 + ^31(^32 -t2l)q2 +0(At3)]. 


If ^32 — ^21 = ts + ti — 2^2 = 0, implying that the interpolation for d? / di? is 
done at the central value t2, then 

s(qi,q3) = ^hih2hiPi x pg •q2 (i + o(At2)) ; 

else, if t2 ^ [h + t3)/2 the last factor is just (1 + 0{At)). Using ^ in ^ 

^(qi,q2,q3) = pi x pa • |-^<2ii32i3iq2+| 



'77^21^32^31(^32 ^ ^21) 
6 



3(q2 • q2)q2 



2 
95 



q2 



If, as above, ^32 + ti2 = tg + ii — 2i2 = then 



^(qi,q2,q3) 

and we can conclude 
To compute P we need 



1*21*32*31/01 X P3 • q2 (1 + ©(At^)) 



-- = l + 0(At2) 



dp 



dt 



(24) 



to make a Taylor expansion of pi in ^2 

P» = P2 + t,2??V2 + -^{-ifp2 + r/V2 + Kr/2n2) + 0(Ai3). 
This implies that 

Pi X p3 • p2 = ^ [*i2f?V2 X tl^mfn2 - t32??V2 X tl^nrj^ na] • P2 + 0{At^) 



and the ©(At ) term vanishes thus 



P = -^{tuth - *32<?2) (1 + 0{At^)) = ^*21*32*31 (1 + 0{At')) 
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and 



Co 



B ^pi xp3.q2(l + 0(Ai))' 
In the denominator pi x pa computed to order At^ is 



If t 



32 



Pl X P3 = 

^21 = ^3 + tl 

Co = 



32 



^31 "2 

- 2t2 = then 
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(?7n2 - KTy^ V2) + C(A<^ 



Aii3i ?yg2q2 • n2 + 0(At3) ^q2-n2 



(1 + (OA<2)) 



(25) 



(26) 



otherwise the last factor is (1 + ©(At)). 

We can conclude that if the topocentric correction is neglected the coeffi- 
cients of the two dynamical equations ([5|) and ((2T|) are the same to zero order 
in At and also to order 1 if the time ^2 is the average timcF^I. 



3.2 Topocentric Correction in Laplace's Method 

Now let us remove the approximation that the observer sits at the center of 
the Earth and introduce the topocentric correction into Laplace's method. The 
center of mass of the Earth is at q0 but the observer is at q = q® + P. Let 
us derive the dynamical equation by also taking into account the acceleration 
contained in the geocentric position of the observer P(t) such that 



Multiplying by -n and using eq. ^ 



ma 



prfn 



q®- 



9©- 



p- 



P n 



P n 



The term PP ■ h/r^ can be neglected. This approximation is legitimate be- 
cause P/q^ < 4.3 X 10~^ and the neglected term is smaller than the planetary 
perturbations. Thus 



C^ = (l-A„)-% 



where 



C 



A„ 



P n 



g©P ■ _ 

Mq©-n (Ai/g2 )qgj . 



(27) 



(28) 



^"This equivalence is taken for granted by many authors but only in [Poincar^ "1906] we 
have found the basic idea of the computations above. 
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Note that A„ is singular only where C is also singular. The analog of eq. ©J 
again neglecting Oip/q^), is 



pf] + 2/577 



4 



( 



) 




(29) 



The important fact is that A„ and A^, arc by no means small. The centripetal 
acceleration of the observer (towards the rotation axis of the Earth) has size 

i?0 cos6' where fi® is the angular velocity of the Earth's rotation, the 
radius of the Earth and 6 the latitude; the maximum of ~ 3.4 cms~^ occurs 
at the equator. The quantity fJ-fq^ in the denominator of A„ is the size of the 
heliocentric acceleration of the Earth, ~ 0.6 cm s^^. Thus |A„| can be > 1, and 
the coefficient 1 — A„ very different from 1; it may even be negative. This leads 
to the conclusion that without taking into account the topocentric correction the 
classical method of Laplace is not a good approximation in the general cas 

The common procedure when using Laplace's method is to apply a negative 
topocentric correction to go back to the geocentric observation case. However, 
in doing this some value of p is assumed as a first approximation. If this value is 
approximately correct, by iterating the cycle (topocentric correction - Laplace's 
determination of p) convergence is achieved. If the starting value is really wrong 
the procedure may well diverge, {e.g. a value p = I AU is assumed by default 
and the object is actually undergoing a close approach to the Earth). Moreover, 
in this way the information contained in the parallax is not exploited in an opti- 
mal way. Unless a better way is found to account for the topocentric correction 
there are reliability problems discouraging the use of Laplace's method when 
processing a large dataset, containing discoveries of objects of different orbital 
classes and therefore spanning a wide range of distances. 

3.3 Gauss-Laplace equivalence, Topocentric 

When taking into account the displacement P the Taylor expansion of qi{t) of 
eq. ([23)1 is not applicable. We need to use 



where q2(t) and its derivatives contain also P{t). By using eq. (p6)) and assuming 
t2i = h2, cq. (HH) and ^ become 



^^When observations from different nights arc taken by the same station at the same sidereal 
time the topocentric correction in acceleration cancels out. In this case the classical Laplace 
method is a good approximation. 




^(qi,q3) = ^ t2lt32tli n2 •q2 +0(At') 
3 

^(qi, q2, qs) = ^ t2lt32tl, n2 • q2 + 0(At6) . 
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Note that q2 does not appear in A at this approximation level. Thus 

A _ g| n2 ■ q2 + OjAt'') 
° B Mn2-q2+0(Ai2) 

and once again neglecting P/qt^ terms 

92 "2 • qe2 ql "2 • P2 



^ n2 • q2 ^ n2 • q2 

^ gj "2 ■ P2 

9e2 M n2 • q2 



Finally 



then 



n2 • q2 = 92 n2 • \ = 92 n2 • qe2 + [ — 

92 92 / V V 92 



ho = l~ '^^^"^•^^ + 0(At2) + o(^]=l- A„2 + OiAf) + O 

/U n2 • q2 V 92 / V 92 

where A„2 is the same quantity as A„ given by cq. (j28p and computed at t = ^2- 
The conclusion is that Gauss' method used with the heliocentric positions 
of the observer q^ is equivalent to Laplace's method with topocentric correction 
to lowest order in At and neglecting the very small term 0{P2/q2)- 

3.4 Problems in Topocentric Laplace's Method 

The results obtained in this Section can be summarized as follows: contrary to 
common belief, Gauss' method is not equivalent to Laplace's unless Gauss' is 
artificially spoiled by not using the observer position in eq. and ([^0]) . In 
its rigorous form, Gauss' method accounts for the topocentric correction with 
a consistent approximation. The question arises whether we could account for 
the topocentric correction in Laplace's method (without iterations) by adding 
the term A„ from eq. (j28p . Surprisingly, the answer is already contained in the 
literature in a 100 year old paper by a famous author [Poincare 19061 pag. 177- 
178]. To summarize the argument of Poincare, plots showing the shape of the 
topocentric corrections as a function of time and a short citation are enough. 

Figure [2] shows the simulated path of an approaching NEC as seen from 
an observing station (in this example in Hawaii). The darker portions of the 
curve indicate possible observations that have an altitude > 15°. The overall 
apparent motion of the asteroid from night to night cannot be approximated 
using parabolic motion segments fitted to a single nighlF^. For the geocentric 
path the parabolic approximation to p(t), used by Laplace, would be applicable. 



'^^Our translation of Poincare: It is necessary to avoid computing these quantities by starting 
from the law of rotation of the Earth. 
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Figure 2: The path in the sky of an approaching NEO. This example is (101955) 
1999 RQ36 as it would have been seen in July 2005 if an observatory on Mauna 
Kea had been observing continously. The curve resembling a parabola gives 
simulated observations from the geocenter. 

Figure [3] shows graphically that topocentric observations contain informa- 
tion beyond what is contained in the average angles and proper motion (the 
attributable, see Section [5|). Thus, to reduce the observations to the geocenter 
by removing the topocentric correction is not a good strategy. 

Poincare suggests computing what we call A„ by using a value of P obtained 
by interpolating the values P(ii) at the times ti of the observations which are 
not limited to 3 (one of the advantages of Laplace's method). We have im- 
plemented Poincarc's suggestion to improve Laplace's method in our software 
system OrbFiQ but we still need to test it properly. 

When the observations are performed from an artificial satellite (such as the 
Space Telescope or, in the future, from Gaia) the acceleration P ~ 900 cm 
and the A„ and coefficients can be up to ~ 1, 500. A few hours of observations 
extending to several orbits can produce multiple kinks as in IMarchi et al. 2004|. 
Figure 1] containing important orbital information. 

4 Qualitative Theory, Topocentric 

In rectangular heliocentric coordinates {x, y) where the x axis is along q2 (from 
the Sun to the observer) we have p2 = v'g| + a;^ + j/^ — 2x?2 and 7'2 = \/ x'^ + y^, 

^^In Version 3.4.2 and later; see |http : //newton.dm.unipl ■ It/orbf lt/| 



14 



15 




15' ' ' ' ' ' ' ' ' 

-40 -30 -20 -10 10 20 30 40 

Right Ascension difference (arcsec) 



Figure 3: The same data as in the previous figure after removing tlie best fitting 
linear functions of time in botli coordinates. In this case the curves represent the 
content of information beyond the attributable. The larger loop is from Mauna 
Kea while the small curl near (0, 0) is for a geocentric observer. Coordinates 
are differences in RA and DEC in radians. 



thus we can consider the function 



(a-2+y2)3/2^ 

The dynamical equation eq. (PT|) can be seen as describing the level lines Co = 
const in a bipolar coordinate system (r2, ^2)- Note that Co = is the zero circle 
r = ro = q/ \/ho for Hq > and is otherwise empty. A simple computation of the 
partial derivatives of Co shows that the only stationary points of Co are the pairs 
(a;, y) with y = and x a solution of the equation (/lol^^P — 92)^ ^ •^Q2i^ ~ ^2). 
For /lo < it has only one solution, xi, with < xi < (72. For ho > there is 
always at least one solution x < —tq < 0. If < /iq < 1 there are two additional 
solutions, xi and X2 such that < xi < q2 < tq < X2- For ho > 1 there are no 
positive solution^. 

The function Co(a;, y) has a pole of order 3 at (a;, y) — (0, 0), with limr^^o C'o — 
—00, and a pole of order 1 at ((72,0) with linXr^^oC'o ~ +00 for > I and 
= —00 for ho < 1. For ho = 1 there is at ((72, 0) a more complicated singularity: 
as shown by Figure [T] there is no unique limit value for p2 0. 

'^■'The quantity C appearing in the topocentric Laplace's method defines exactly the same 
function of (x, y), with hg replaced by 1 — A„. 
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Level curves of with h^<0 




-0.5 0.5 1 1.5 

Figure 4: Level curves of Co{x,y) for ho = —0.5. Note there is no zero circle. 

4.1 Topology of the level curves of Co{x,y) 

The qualitative behavior of the level lines of Co{x,y) is different in the three 
cases ho <0 (Figure |4|), < /iq < 1 (Figure E]) and ho > 1 (Figure [H). 

The number of solutions of the dynamical equation {i.e. along a fixed 
topocentric direction) can be understood by evaluating the degree 8 polyno- 
mial (|22|1 on the zero circle 

^(-o) = C^o -Ca (l - '^o^') • 



Table 1: The number of preliminary orbit solutions, the columns give: [1] the 
number of preliminary orbit solutions (for a given C and e), [2] the number of 
positive roots of the polynomial equation (p2)) . [3] the number of spurious roots. 
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1 or 3 
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C < 
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1 or 3 
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C > 




1 or 3 




1 or 3 


or 2 



We summarize the possible numbers of solutions, for a given direction of 
observation e, in the different cases, depending upon the value of Hq and the 
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Figure 5: Level curves of Co{x, y) for ho = 0.5, including the zero circle (dashed). 

sign of C, in Table [TJ Note that in the Table we are not making the assumption 
of Charlier, that some solutions must exist, for the reasons given in Section [2721 
By spurious we mean a root of the polynomial (|22p corresponding to p < in eq. 
(Plj). This is not a complete qualitative theory replacing Charlier's for /iq = 1, 
but already shows that the number of solutions can be quite different from the 
classical case e.g. , 2 solutions near opposition and up to 3 at low elongation. 
For a fully generalized qualitative theory see [Gronchi 2007] . 

4.2 Examples 

We would like to find examples in which the additional solutions with respect to 
the classical theory by Charlier are useful. That is, cases in which the additional 
solutions provide a preliminary orbit closer to the true orbit and therefore more 
suitable as a first guess for the differential corrections procedure. 

An example in which there arc two solutions while observing in a direction 
close to the opposition is shown in Figure [T] Of the two intersections of the 
observing half line with the relevant level curve of Co, the one leading to a 
useful preliminary orbit is the nearer one which has p2 = as counterpart in 
the classical theory. The farther one leads to a preliminary orbit with e ~ 10. 
We have used the formulae of [Milani et al. 200"4] to compute the maximum 
possible p2 along the observation direction compatible with e < 1. 

Another interesting feature of this example is that the preliminary orbit 
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Level curves of with h^>^ 
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Figure 6: Level curves of Co (a;, y) for ho = 1.1, including the zero circle (dashed). 

using the nearer solution has residuals of the 6 observations with RMS = 66 
arcsec while the one using the farther solution has RMS = 2.5 arcsec. This 
implies that if only one preliminary solution were passed to the next processing 
step by selecting the one with lowest RMS the good solution would be discarded. 

To find a significant example with 3 solutions is not easy because in many 
cases the third solution, the nearest to the observer, has p2 too small for the 
heliocentric 2-body approximation to be applicable. A value p2 < 0.01 AU 
corresponds to the sphere of influence of the Earth, i.e., the region where the 
"perturbation" from the Earth is actually more important that the attraction 
from the Sun. Thus, a solution with such a small p2 must be discarded, because 
the approximation used in Gauss' and Laplace's method is not valid. 

To show how our arguments on the number of solutions applies to a real 
case (as opposed to a simulation as in the example above) we have selected 
the asteroid 2002 AA29 and used observations from the first three nights (9, 11 
and 12 January 2002). With the values Cq ~ 1.653 , ho = 1.025 we obtain 
from the observations and an elongation ~ 111° there is only one solution with 
P2 — 0.0045 (see Figure[51 left), which easily leads to a full least squares solution 
with p2 = 0.044. Although the value of ho is not very far from 1 the existence 
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Figure 7: An example with two solutions near opposition: for Hq ~ 0.613 
the direction of observation (dotted) has two intersections with the level curve 
Co{x,y) = 0.4 (continuous); the zero circle is dashed. The positions in the 
observation direction with a bounded orbit are drawn as a continuous line. 



of the solution depends critically on ft-o ^ 1 7^ 0. If the value of ho had been set 
to 1 we would find no solution (see Figure [8l right) . 

4.3 Implementation issues 

We need to implement the algorithms discussed in this paper for the computa- 
tion of preliminary orbits in a way which is suitable for a large observation data 
set; we need to satisfy three requirements. 

The first requirement is to obtain the solutions to the polynomial equations 
such as (j22p in a way which is fast and reliable in providing the number of 
distinct real solutions. In this way we can fully exploit the understanding on 
the number of solutions (with topocentric correction) which we have achieved 
in this section. This is made possible by the algorithms computing the set of 
roots of a polynomial equation at once (as a complex vector) and with rigorous 
upper bounds for the errors including the ones generated by roundoff. We use 
the algorithm by |Bini 1996] and the corresponding public domain software!^. 

The second requirement is to improve the preliminary orbit as obtained 
from the solutions of the degree 8 polynomial equations in such a way that 
it is as close as possible to the "true" solution to be later obtained by differ- 
ential corrections. There is such an immense literature on this topic that in 
this paper it is not even appropriate to discuss the references. Conceptually, 
as shown by [Celletti and Pinzari 2005] , each step in the iterative procedures 
used in differential corrections can be shown to increase the order in At of the 

^^For the Fortran 77 version |http: //wot. netllb.org/mimeralgo/nalO[ For Fortran 90 
|http : //users .blgpond. net ■ au/amiller/pzeros . f 90| ■ 
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Figure 8: For the preliminary orbit of 2002 AA29 the relevant level curve (Co = 
1.653) is shown (continuous) in the same plane of Figure [1] the zero circle 
(dashed) and the observation direction (dotted) are also shown. Left: using the 
actual value ho = 1.025. Right: using a value of /iq = 1 that does not account 
for the topocentric correction. 

approximation to the exact solutions of the 2-body equations of motion. How- 
ever, jCelletti and Pinzari 2006] have also shown that an iterative Gauss map 
can diverge when the solution of the degree 8 equation is far from the fixed 
point of the iterative procedure, outside of its convergence domain. 

We have implemented one of the available iterative improvement algorithms 
for Gauss' method and have found that it provides in most cases a preliminary 
orbit much closer to the least squares solution and therefore a more reliable first 
guess for the least squares algorithms. We have also found that the Gauss map 
diverges in a small fraction of the test cases but still often enough to significantly 
decrease the efficiencj0 of the algorithm. In some cases the number of orbits 
for which the Gauss map converges is less than the number of solutions of the 
degree 8 equations. It can happen that one of the lost degree 8 solutions was 
the one closest to the "true" orbit and the only one which can be used to obtain 
the best least squares solutions. One method to obtain the highest efficiency 
without an inordinate increase in the computational cost is to run two iterations, 
one with and one without the Gauss map. 

The third requirement is to use modified differential corrections iterative 
algorithms with larger convergence domains in such a way that even when the 
geodetic curvature (contained in the coefficients C and Co of the two methods) 
is poorly constrained by the available observations (because the arc length on 
the celestial sphere is too short) the very rough preliminary orbit solution can 
lead to a least squares solution. This is discussed in the next section. 

^^See Section |6] for the definition of the metric. 
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5 Weak preliminary orbits 



An essential difference between the classical works on preliminary orbits and the 
modern approach to the same problem is that the effects of the astrometric errors 
cannot be neglected, especially in the operating condition of modern surveys: 
they use shorter observed arcs, thus the deviations of the observed path on the 
celestial sphere from a great circle may not be significant. 



5.1 Uncertainty of Curvature 

The explicit computation of the two components of curvature of interest for orbit 
determination, geodesic curvature k and along track acceleration can be per- 
formed by using the properties of the orthonormal frame ([T]) by straightforward 
computation using the Riemannian structure of the unit sphere [Milani et al. 2006b|, 
Section 6.4]. The results are 



I (J d — d (5) 



cos 5 + a 



+ {6f sin = K(a, 5, d, j, d, 5) (30) 



d d cos^ 5 + 55— (d)^ 5 cos 5 sin 5 — ^, d, j, d, 5) . (31) 



Given these explicit formulae it is possible to compute the covariance matrix of 
the quantities (k, t)) by propagation of the covariance matrix of the angles and 
their derivatives with the matrix of partial derivatives for k and f] 



9(k,?7) 



d{a, 5, d, 5, d, 5) 



d{K,ri) 



1 T 



d{a, 5, d, (5, d, 5) 



(32) 



The covariance matrix T^^s for the angles and their first and second derivatives 
is obtained by the procedure of least squares fit of the individual observations 
to a quadratic function of time. The partials of k and f] are given below (note 
that the partials with respect to a are zero). 
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-6 a + d? cos 5 sin (5 +a5 
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a cos (5 



di) 

as 



Note that the last four of the partials above, the 2x2 matrix d{K,r])/d{d, S), 
contribute to the principal part of the covariance of (k, ?)) for short arcs as 
discussed in the next subsection. 

We use a full computation of the covariance matrix without approximations 
to assess the significance of curvature by using the formula from [Milani et al. 2006b| 



(33) 
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and we assume that the curvature is significant if > X 



2 

min 



5.2 The Infinite Distance Limit 

The problem of low values of C can occur in two ways: near the zero circle 
and for large values of both p and r. On the other hand, the uncertainty in 
the estimates of the deviations from a great circle will depend upon the length 
of the observed arc (both in time At and in arc length ^ '^Ai). For short 
observed arcs it may be the case that the curvature is not significant. Then the 
preliminary orbit algorithms will yield inaccurate preliminary orbits which may 
fail as starting guesses for differential corrections. 

We will now focus on the case of distant objects. We would like to estimate 
the order of magnitude of the uncertainty in the computed orbit with respect to 
the small parameters r, b where u is the astrometric accuracy of the individual 
observations (in radians) and t = n^At, b ~ ^ are small for short observed 
arcs and for distant objects, respectively. Note that the proper motion rj for 
6^0 has principal part n^b - the effect of the motion of the Earth. The 
uncertainty in the angles (a, 5) and their derivatives can be estimated as follows 

The uncertainty of the curvature components (k, -ff) should be estimated by the 
propagation formula ((32|) but it can be shown that the uncertainty of ((5, d, 5) 
contributes with lower order terms. Thus we use the estimates 

d{a,5) 

and obtain 

_ r o(&-v-2) o(&-v-2)n| 1 



0{b-^)n-' 0{b-^)n-' 
Oil) 0(1) 
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To propagate the covariance to the variables {p, p) we use the equation, ob- 
tained by ehminating r between eq. (O and (|2ip . an impHcit equation connecting 
C and p 



F{C,p) = C ^ 
we have C b^^ 



ill 



2(70pcose)'^/2 



1 + A, 







(34) 



For 6 — > we have C 1 and thus C ^ and is of the same order as the 

small parameter b. Although C depends upon all the variables (a, S, d, S, d, 6), 
its uncertainty mostly depends upon the uncertainty of k and thus, ultimately, 
upon the difficulty in estimating the second derivatives of the angles. Next, we 
compute the dependence of the uncertainty of (p, p) upon the uncertainty of 
(k, ?)). From the derivatives of the implicit function p{k), assuming cos e, 77, n to 
be constant and keeping only the term of lowest order in q/ p, 
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In the same way from ((29)) we deduce 77 
partial derivatives 

1^ = Tie % 0{1) , 
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For the covariance matrix 
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we compute the main terms of highest order in 6 ^ , t ^ as 
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-3 ^-2 
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(35) 



In conclusion, if the variables (p, p) are measured in the appropriate units (AU 
for p and AU for p) the uncertainties are of the same order and there is no 
reason to suppose that one of the two will be better determined than the other. 

This conclusion is different from the one of [Bernstein and Khushalani 2000] . 
We are making no assumption about the orbit, just on the distance of the ob- 
served object, e.g. the proper motion may not be well aligned with the ecliptic 
as is the case for a low eccentricity, low inclination "classical TNO" . This is due 
to the fact that our main concern is reliability. We do not want to use an orbit 
computation method which might preferentially fail on unusual orbits, e.g. for 
a long period (even hyperbolic) comet discovered at large distance. Wc do agree 
with [Bernstein and Khushalani 2000] on the fact that for a TNO observed only 
over an arc shorter than one month there is very often an approximate degener- 
acy that forces the use of a constrained orbit (with only 5 free parameters). We 
only claim that the weak direction, which is essentially in the (p, p) plane, may 
vary and is generally not along the p axis [Milani et al. 2005b[ Figures 3-6] . In 
Section [6. II we confirm this argument by a numerical test on TNOs. 
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5.3 Prom Preliminary Orbits to Least Square Solutions 

The procedure to compute an orbit given an observed arc with > 3 nights of data 
(beUeved to belong to the same object) begins with the solution of the degree 
8 equation ([2^ and ends with the differential corrections iterations to achieve 
a full least squares orbit (with 6 solved parameters). However, for algorithms 
more efficient than the classical ones there are up to four intermediate steps. 

We need the definition of Attributable [Milani et al. 200T| : the set of 4 vari- 
ables (a, S, a, S) estimated at some reference time, e.g. ^2, by a fit to the obser- 
vations. It is possible to complete an attributable to a set of orbital elements by 
adding the values of range and range rate {p2, P2) at the same timj^. For each 
attributable we can determine an Admissible Region which is a compact set in 
the i^piiPi) plane compatible with Solar System orbits [Milani et al. 2004] . 

The optional intermediate steps are 

1. an iterative Gauss map to improve the solution of the degree 8 equation; 

2. adding to the preliminary orbit (s) another one obtained from the At- 
tributable and a value for {p2,P2) selected inside the Admissible region 
(see details below); 

3. a fit of the available observations to a 4-parameter attributable at time ^2; 
the values of p2 and p2 are kept fixed at the previous values; 

4. a fit of the available observations constrained to the Line Of Variations 
(LOV), a smooth line defined by minimization on hyperplanes orthogonal 
to the weak direction of the normal matrix. 

Intermediate step 1 has been discussed in Section l473l 

For intermediate step 2 we need to distinguish two cases depending upon the 
topology of the Admissible Region (Milani et al. 2004] . If it has two connected 
components (this occurs for distant objects observed near opposition) we select 
the point which is the center of symmetry of the connected component farthest 
from the observer. This corresponds to an orbit with e < 1 although it is not, 
in general; a circular orbit which may be incompatible with the Attributable. 

If the Admissible Region is connected then we select the point along the 
symmetry line p2 ~ const at 0.8 times the maximum distance p2 compatible 
with e < 1- This case always occurs near quadrature; if the object is indeed 
distant, thus has a low proper motion 77, the selected point is also far. 

In any case, the selected point (p2,P2) in the admissible region completed 
with the Attributable provides an orbit which is compatible with the given 
Attributable and belongs to the Solar System; this is called a Virtual Asteroid 
VA) [Milani 2005] . The VA method provides an additional preliminary orbit. 
This does not matter when there are already good preliminary orbits computed 
with Gauss' method (with significant curvature). We shall see in Section [Ol 

^'^The epoch time of the elements is to = ^2 — P2/c; these elements are said to bo in 
Attributable Coordinates [Milani et al. 2005b) . 
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that for TNOs this additional prehminary orbit is required in many cases, the 
majority of cases near quadrature, because the curvature is hardly significant. 

Intermediate Step 3 is essentially the method proposed by D. Tholen and also 
available in his public domain software KNOBS. It has already been tested in 
the context of a simulation of a next generation survey in [Milani et al. 2006a| . 

Intermediate step 4 is described in full in [Milani et al. 200"5a| . Our preferred 
options are to use either Cartesian Coordinates or Attributable Elements scaled 
as described in [Milani et al. 2005al Table 1]. 

The steps listed above are all optional and indeed it is possible to compute 
good orbits in many cases without some of them. However, if the goal is a very 
reliable algorithm it is necessary to use them with a smart connecting logic. As 
an example. Step 1 is used in a first iteration, can be omitted in a second one. 
Step 2 is essential for distant objects. Step 3 is used whenever the curvature 
is not significant i.e. when the observed arc is of type 1 [Milani et al. 2006bj 
which can be tested e.g. by eq. (|33l) . Step 4 is important for weakly determined 
orbits, otherwise the differential corrections may diverge when starting from an 
initial guess with comparatively large residuals. Even then. Step 4 may fail 
and, in turn, diverge under differential corrections. In this case the differential 
corrections restart from the outcome of the previous step. This connecting logic 
is an extension of the one presented in [Milani et al. 200"5al Figure 5] . 

6 Tests 

The tests we are using are obtained by running a simulation based upon a 
Solar System Model - a catalog of orbits of synthetic objects as described in 
[Milani et al. 2006a| . Given an assumed observation scheduling and instrument 
performance we compute the detections of the catalog objects above a threshold 
signal to noise ratio and record the corresponding simulated astrometric obser- 
vation including astrometric error. In the simulations used here we have not 
included false detections (not corresponding to any synthetic object). 

Then we assemble into tracklets the detections (from the same observing 
night) which could belong to the same object. The tracklets are assembled in 
tracks from several distinct nights (in this context, at least three nights are re- 
quired) . For these simulations we have used the algorithms of [Kubica et al. 2007) 
to assemble both tracklets and tracks. 

When the number density of detections per unit area is low both tracklets 
and tracks are (almost always) true i. e. they contain only detections of one and 
the same synthetic object. When the number density is large, as expected from 
the next generation surveys, both tracklets and tracks can be false (containing 
detections belonging to different objects). This is why a track needs to be 
confirmed by computing an orbit: first a preliminary orbit, then by differential 
corrections another orbit which fits all the observations in the least squares sense. 
The structure containing the track and the derived orbit with the accessory data 
necessary for quality control (covariance matrix, weights, residuals, statistical 
tests) is called an identification [Milani et al. 200"7) . 
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The purpose of the tests is to measure the performance of the algorithms 
described in this paper, according to the following criteria: 

• Efhcicncy E: the fraction of true tracks for which good preliminary/least 
squares orbits were calculated. 

• Accuracy A: the fraction of returned orbits that correspond to true tracks. 
I.e., the orbit computation should fail on false tracks (cither no preliminary 
orbit or no least squares orbit or residuals too large) . 

• Speed S: Reciprocal of the CPU processing time. 

• Goodness G: the fraction of least squares orbit close enough to the ground 
truth orbit of the object to allow later recovery {e.g. in another lunation). 

The Speed criterion is less important than the others for the reasons ex- 
plained in Section [TJ Nevertheless, we need to check that the very large data 
sets expected from the next generation surveys can be processed with reasonable 
computational resources. 

Note that these tests should not be confused with tests of the performance 
of next generation surveys like Pan-STARRS or LSST. The purpose is to show 
that whatever the rate at which new objects arc observed their discoveries will 
not be lost because of inefficiency in the orbit determination procedure. 

6.1 Small targeted test 

Since the orbits of MBA and Jupiter Trojans are easier to compute than NEOs 
and more distant objects [Milani et al. 200"6a] . to assess the Efficiency of these 
algorithms we have prepared four targeted simulations: two containing only 
observations of NEOs and two with TNOs only. In both cases, one simulation 
uses a surveying region near opposition and the other surveys the so called sweet 
spots at solar elongations between 60° and 90°. The metrics being measured by 
these tests are EfRciency and Goodness: Speed is irrelevant for such small data 
sets and Accuracy is not a serious issue because the number density per unit 
area is small (indeed. Accuracy is 100% in all tests of this Subsection). 

The first part of Table prefers to the simulation including only NEO around 
opposition. Note that the Incomplete Identifications are due to the fact that 
the algorithm used to assemble the tracks operates on observations belonging to 
the same lunatioiF^ while the simulation included two lunations. The separate 
orbits obtained in the two lunations for the same object can be joined later with 
other algorithms, e.g., the ones of [Milani et al. 200T| . Thus there are only 4 
"failures" , true tracks which have not been confirmed by the orbit computation. 
The lower part of Table [2] refers to the simulation including only NEOs at the 
sweet spots. There are only 8 3-nighters without an orbit. 

In conclusion, in both NEO simulations the Efficiency is very high but not 
perfect, especially at the sweet spots. Note that in such small tests it is always 

^*An arc of one month has an excessive curvature jKubica et al. 20071 Section 8.2]. 
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Table 2: For each of the NEO simulations, separately for objects observed on 
a different number of nights, the columns give: [1] Total number of objects, [2] 
Number of Complete Identifications (containing all the tracklets belonging to 
the object), [3] Efficiency (defined as [2] /[I], in %), [4] Number of incomplete 
Identifications, [5] Fraction [4] /[I] in % of incomplete Identifications, [6] Number 
of objects lost (no confirmed Identification), [7] Fraction [6]/[l] in % lost. 





[1] 


[2] 


[3] 


[4] 


[5] 


[6] 


[7] 


Observed 










Inc. 




Lost 


Arc 


Total 


Compl. 


Effic. 


Inc. 


Fraction 


Lost 


Fraction 


Opposition 
















3-nighters 


1123 


1119 


99.6% 





0.0% 


4 


0.4% 


4-nighters 


2 


1 


50.0% 


1 


50.0% 





0.0% 


6-nighters 


123 





0.0% 


123 


100.0% 





0.0% 


Sweet Spots 
















3-nighters 


397 


389 


98.0% 





0.0% 


8 


2.0% 


6-nighters 


63 





0.0% 


63 


100.0% 





0.0% 



possible to increase the Efficiency to 100% by running additional iterations 
with increased computational intensity. However, this would not provide useful 
indications on what should be done with a much larger data set (e.g., we could 
increase Efficiency at the expense of Accuracy). Nevertheless, it is useful to 
examine the few failure cases in order to learn about either the limitations of 
our theory or defects in our implementation. None of the cases in which an 
orbit was not computed, although a true track was proposed, result from the 
failure of the differential corrections. They all resulted from a failure of the 
preliminary orbit in the sense discussed in Section 15. 3[ in most cases because 
the degree 8 equation has only invalid roots (either spurious, i.e. p < 0, or p 
positive but small, resulting in a poor 2-body approximation); the VA method 
was not of any help, as expected since it is intended for low curvature cases. In 
other cases some preliminary orbit could be found but the RMS of the fit was 
comparatively large, in the range between 200 and 300 arcsec. 

The upper part of Table [3] refers to the simulation with TNOs around opposi- 
tion. The Incomplete identifications for 3-nighters are cases in which more than 
one tracklet was available on some night. For the cases for > 6 nights the track 
was not proposed for the same reasons discussed in the NEO case above. Thus 
the orbit determination has failed only in one case. Moreover, this single failure 
is not due to the preliminary orbits. The differential corrections stage computed 
a nominal least square orbit that was then refused by the quality control stage 
not because of small residuals (RMS 0.084 arcsec) but due to systematic trends 
(e.g., a slope) with a signal to noise ~ 2.5. 

The lower part of Table [3] refers to the simulation including only TNOs in 
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Table 3: For each of the TNO simulations the same data as in Table [2] 





[1] 


[2] 


[3] 


[4] 


[5] 


[6] 


[7] 


Observed 










Inc. 




Lost 


Arc 


Total 


Compl. 


Effic. 


Inc. 


Fraction 


Lost 


Fraction 


Opposition 
















3-nighters 


2005 


2001 


99.8% 


3 


0.15% 


1 


0.05% 


4-nighters 


18 


18 


100.0% 





0.00% 





0.00% 


5-nighters 


3 


3 


100.0% 





0.00% 





0.00% 


6-nighters 


670 





0.0% 


670 


100.0% 





0.00% 


7-nighters 


13 





0.0% 


13 


100.0% 





0.00% 


8-nighters 


1 





0.0% 


1 


100.0% 





0.00% 


9-nighters 


6 





0.0% 


6 


100.0% 





0.00% 


Sweet Spots 
















3-nighters 


2493 


2491 


99.9% 





0.00% 


2 


0.08% 



the sweet spot regions. Here the simulated scheduling included only 3 nights 
and there are only two cases in which an orbit was not computed. As in the 
opposition simulations these failures are due to tight quality control thresholds 
because of systematic trends with signal to noise between 3 and 5.5. 



Table 4: Fraction of objects lost 3-nighters using different algorithms. 



Simulation 


Best 


1 Pre 


1st It. 


No VA 


No 4fit 


No LOV 


LSQ 


NEO 0pp. 


0.40% 


0.50% 


3.2% 


0.4% 


0.4% 


0.4% 


0.4% 


NEO Sw. 


2.00% 


13.4% 


18.1% 


2.3% 


2.3% 


2.8% 


2.8% 


TNO 0pp. 


0.05% 


0.05% 


0.05% 


38.3% 


38.4% 


45.7% 


47.7% 


TNO Sw. 


0.10% 


0.10% 


0.10% 


75.3% 


75.3% 


82.6% 


82.8% 



In conclusion, the preliminary orbit algorithms have not shown even one case 
of failure in the TNO simulations. However, we need to assess the proportion of 
this success due to the improved but classical method of Gauss rather than to the 
Virtual Asteroid method. The latter is expected to be especially effective for the 
low curvatures typical of TNOs. We have thus run again the simulations with 
a version of the software not containing the VA method. The difference of the 
results with the previous ones measures the contribution from the VA method as 
shown in Table [4] in the column labeled "No VA" . In the opposition simulations 
giving up the VA method resulted in > 1/3 of the 3-night TNOs being lost 
due to a lack of orbit computation. In the sweet spot simulations the same 
situation occurred for 3/4 of the 3-night TNOs. The reason for this difference 
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is easily understood as near quadrature the TNO have an even smaller proper 
motion than at opposition and the curvature is very often not well measured. 
The conclusion is that the VA method is essential for TNOs while it is almost 
irrelevant for NEO (only a small contribution in the sweet spots). 

For NEOs the most relevant question is the utility of all the care we have 
exercised in ensuring that no identification is lost because of double (or even 
triple) solutions of the preliminary orbit equations. For this we have run a simu- 
lation in which only one preliminary orbit was passed to differential corrections 
independent from the number of solutions to the degree 8 equation. The selected 
preliminary orbit was the one with lowest RMS of residuals. The results are in 
Table m in the column "1 Pre" which clearly show that to pass to differential 
corrections double (or possibly triple) solutions near quadrature is essential for 
top Efficiency for NEOs. At opposition it matters only in rare caseJ^. 

Another test has been to stop after the first of the two iterations (see Sec- 
tion [531), the one with a tighter control in the RMS of the residuals for the 2- 
body preliminary orbit (set at 10 arcsec in these tests) and using the Gauss map. 
From the results (column "1st It.") it is clear that the second iteration (with 
RMS of preliminary orbit up to 100 arcsec and the solution of degree 8 equation 
directly passed to differential corrections) has no effect at all on TNOs but is 
relevant for NEOs especially in the sweet spots. This is because the NEOs ob- 
served over 3 well spaced nights are observed arcs of type 3 jMilani et al. 2006b] . 
i.e. the information contained in the observations is much more than just (/c, r)) 
of the entire arc which is used in the preliminary orbit. On the other hand, the 
planetary perturbations neglected in the preliminary orbit algorithms are large 
with respect to the observation accuracy (assumed to be 0.1 arcsec). 

Although this paper is mostly about preliminary orbit algorithms we need to 
also assess how much the improved differential corrections algorithms (discussed 
in Section 15. 3[) have contributed to the overall success of these simulations. 
Column "No 4fit" reports the results when the 4-parameter fit step was not 
used. The results are essentially identical to the ones of the "No VA" case 
indicating that the two algorithms must be used together for TNOs. 

The column "No LOV" indicates that the step with the 5-parameter least 
squares fit to obtain a LOV solution has a very large effect for TNO: almost 
1/2 of the identifications in the Opposition case and > 4/5 in the Sweet Spots 
case would not be confirmed by a least squares orbit without the LOV solution. 
Indeed, the LOV solution is the only one available for 53.8% of the TNOs near 
opposition and 81.6% of the TNOs in the sweet spots (as opposed to only 0.3% 
of NEO at sweet spot and 0.1% at opposition). This is not a surprise. Indeed, 
the 3-night TNOs arc almost always observed arcs of type 2 (in some cases 
even type 1) which, according to the tests in |Milani et al. 2006bj . are generally 
not suitable to obtain a well determined orbit. We have also used this test to 
confirm our statement that the weak direction of the LOV is not aligned in one 
special direction. In the opposition test the LOV has an angle (computed with 

^^The first example of Section 14.21 is the only case in which a NEO at opposition would be 
lost by using only the preliminary orbit with lowest RMS as shown by the small change in 
column "1 Pre". 
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the proper scaling as indicated by eq. ([55]) ) between —31° and +17° with the p 
axis. In the sweet spots test it has an angle between —54° and +36° with the p 
axis. In conclusion the weak direction depends upon the elongation. It is closer 
to one of the two axes but a different one in each of the two casej^. 

We get the worst results (column 'LSQ') if neither the 4-fit algorithm or 
the LOV solutions are used and the output of the preliminary orbit is passed 
directly to a full 6-paramcters differential corrections. The difference between 
this column and the one labeled "best" measures the progress between using 
only the classical algorithms and the best combination of algorithms we have 
found so far: for TNOs the difference is very important, with all the steps 
discussed in Section [5.31 essential to achieving the best results. 

6.2 Large scale tests 

The main purpose of a large scale test is to measure the Accuracy. Speed is not 
the limiting factor; besides, the orbit determinations algorithms can be easily 
adapted for parallel processing. Efficiency is not a problem for the overwhelming 
majority of objects, which arc MBA. Accuracy can affect Efficcncy: when there 
are Discordant identifications (with some tracklets in common) if they cannot 
be merged (in an identification with all the tracklets of both) there is no way 
to choose which of the two is the true one. To keep the Accuracy high we then 
have to discard both, which means losing true identifications and decreasing 
Efficiency. On the other hand, we cannot afford low Accuracy since each false 
identification introduces permanent damage to the quality of the resultj^ 

We have prepared simulations for one lunation of a next generation survey 
both near opposition and the sweet spots. The limiting magnitude was assumed 
to be V=24 and the Solar System model was used at full density including 
the overwhelming majority of MBAs and Trojans. Table [5] gives the size of 
the dataset. The focus of this paper is on the objects for which tracklets are 
available in three nights. Objects observed for a smaller number of nights are 
part of the problem in that their tracklets can be incorrectly identified: for such 
large number densities (per unit area on the celestial sphere) false identifications 
happen easily [Milani et al. 2006al figure 3]. 

The first problem concerning Accuracy occurs at the tracklct composition 
stage: some tracklets are false, that is, they mix detections belonging to different 
objects. The question is whether they are identified, thus decreasing Accuracy. 

The second Accuracy problem occurs at the track composition stage. A 
track is just a hypothesis of identification to be checked by computing an orbit: 
at a high tracklet number density most of the tracks are false. The Overhead 
(column marked Overh.) is the ratio between the total number of proposed 
tracks and true ones: it was large, at the sweet spots even above what was 
found in previous simulations [Kubica et al. 20071 Table 3]. 

jBernstein and Khushalani 2000] warn that their arguments arc not apphcable exactly at 
opposition and this is confirmed by our numerical results near opposition. 

False facts are highly injurious to the progress of science, for they often endure long..., 
C. Darwin, The Origin of Man, 1871. 
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Table 5: Simulated data sets: [1] survey region, [2] number of tracklets, [3] 
of which false, [4] number of simulated objects with observed tracklets, [5] of 
which with tracklets in 3 different nights, [6] overhead (see text), [7] objects 
with tracklets in 2 different nights, [8] with tracklets in only 1 night. 



[1] 


[2] 


[3] 


[4] 


[5] 


[6] 


[7] 


[8] 


Oppos. 


654315 


26006 


253289 


164333 


222.8 


41244 


47712 


Sweet sp. 


695067 


59253 


283831 


144903 


501.3 


62177 


76751 



The question is whether the orbit determination stage can produce the true 
orbits with good Efficiency and still reject almost all the false tracks. To achieve 
this, the residuals of the best fit orbits need to be submitted to a rigorous 
statistical quality control^. Our residuals quality control algorithm uses the 
following 10 metrics (control values in square brackets) 

• RMS of astrometric residuals divided by the assumed RMS of the obser- 
vation errors (=0.1 arcsec in these simulations) [1.0] 

• RMS of photometric residuals in magnitudes [0.5] 

• bias of the residuals in RA and in DEC [1.5] 

• first derivative of the residuals in RA and in DEC [1.5] 

• second derivative of the residuals in RA and in DEC [1.5] 

• third derivative of the residuals in RA and in DEC [1.5] 

To compute the bias and derivatives of the residuals we fit them to a polynomial 
of degree 3 and divide the coefficients by their standard deviation as obtained 
from the covariance matrix of the fill^. 

Table 6: Accuracy Results: before and after normalization, total number of 
false identifications accepted, percentage (with respect to the total number of 
identifications), number of identifications containing false tracklets. 

Region All Identifications Normalized 

False % F.Tr. False % F.Tr. 

Oppos. 7093 431 4 80 05 F 
Sweet sp. 1869 1.30 10 29 0.02 



This quality control is currently done with the intervention of a human eye looking at 
the residuals. However, given the number of orbits, the quality control needs to be fully 
automatized. 

^•^When these algorithms are used on real data additional metrics should take into account 
the outcome of outlier removal [Carpino et al., 2003| . For simulations this does not apply. 
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The results are summarized in Tables IH] and [71 As expected, the main prob- 
lem is in Accuracy. Notwithstanding the tight statistical quality controls on 
residuals, while processing tens of millions of proposed tracks a few thousands 
false tracks are found to fit well all their tracklets (columns marked False). The 
numbers are very small with respect to the total number of tracks but they are 
not negligible with respect to the number of true tracks (columns marked %). 
This happens by combining tracklets from 2 (or 3) distinct simulated objects. 
A much smaller number of false tracks contain some false tracklets (columns 
marked F.Tr.). Thus, even the presence of a significant fraction of false track- 
lets does affect neither Efficiency nor Accuracy. 

Given cases like these, with a fit passing all the quality controls, we cannot 
a priori discard any of them since only by consulting the ground truth can we 
know they are false. By further tightening the quality control parameters we 
may remove many false true but some true identifications as well. The values of 
the controls used are already the result of adjustment suggested by experiments 
to find a good trade off between Accuracy and Efficiency. 

The most effective method to remove false tracks is obtained by not consid- 
ering each identification by itself but globally. We have previously defined the 
normalization of lists of identifications in |Milani et al. 2005b[ Section 7] and 
[Milani et al. 2006a[ Section 6]. It removes duplications and inferior identifi- 
cations but also rejects all the Discordant identifications. This is not because 
they are all presumed false, indeed very often one true and one false identifica- 
tion are Discordant, but we do not know which is which. Thus, according to 
our philosophy giving paramount importance to the reliability of the results we 
remove both and sacrifice Efficiency for Accuracy. We have refined the normal- 
ization procedure by checking, for Discordant identification, whether one of the 
two is significantly superior to the other, by comparing the normalized RMS of 
the astrometric residuals: if the difference is more than 0.25 we keep the best. 
The results of this normalization procedure are shown on the right hand side of 
table [H It is clear that the number of false tracks can be reduced to negligible 
values in this way. 



Table 7: Efficiency Results 
Obj.Type Opposition Sweet Spots 

Total Eff.% Eff.No.% Total Eff.% Eff.No.% 



All 


161146 


97.3 


95.9 


144903 


98.0 


97.4 


MBA 


154700 


97.3 


95.8 


135911 


98.0 


97.4 


NEO 


353 


90.4 


90.4 


271 


80.1 


80.1 


Tro 








6894 


97.9 


97.8 


Com 


665 


98.6 


97.6 


253 


98.0 


97.6 


TNO 


5428 


97.7 


97.7 


1574 


98.7 


98.7 



However, the Efficiency also changes as a result of the normalization. In 
Table [7] the rows indicate how the results differ according to the orbital class 
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of the simulated objects. The row marked Com includes Centaurs, long period 
and short period comets. Note that the sweet spots simulation does not detect 
any Jupiter Trojans because the Trojan swarms were not in these directions at 
the time of the synthetic observations. The columns marked Eff.% show the 
Efficiency before and after Eff.No.% normalization. As a rule of thumb, on 
average a little more than 1 true identification is lost for each false rejected. 

Tabic [7] refers to the objects observed with 3 tracklets in 3 nights. A minor 
additional problem occurs in the opposition simulation. For the 3, 187 objects 
observed with > 3 tracklets in 3 nights the proposed tracks may be true but 
incomplete, e.g., have only 3 tracklets when the corresponding object has more. 
In fact, for 13.5% of these objects one or more tracklets fail to be included in 
the identificatiorF^. Although their number is not large these incomplete iden- 
tifications are difficult to be fixed at some later processing stage. The solution 
is to run an algorithm of attribution |Milani et al. 2001] to add the omitted 
tracklets. We have tested this method and found it to be 99.8% efficient. 

The question then is: are the results of Tables [H] and [7] satisfactory and, if 
not, what else can be done? Note that the Efficiency for NEOs and TNOs is 
not affected by normalization (because the number densities arc much lower). 
It might be argued that loosing a few percent of MB As is not important. Never- 
theless, we claim that even this problem can be solved together with the other, 
possibly more important, of tens of NEOs and comets lost for other reasons. 

By separately analyzing the Efficiency of the three steps of the procedure 
(track composition, orbit computation, normalization) we have established that 
the algorithm to generate tracks has been 97.6% efficient at opposition and 
98.7% at the sweet spots. The efficiency of the orbit computation procedure on 
the proposed true tracks has been 99.8% efficient at opposition and 99.3% at the 
sweet spots. The normalization procedure has been 98.6% efficient at opposition 
and 99.4% at the sweet spots. It follows that the different steps are well balanced 
in their performances and there is not much room for improvement. Thus the 
solution is to use a two iteration procedure. 

The normalization procedure generates two outputs: the new list of iden- 
tifications and the list of leftover tracklets which have not been used in the 
confirmed identifications. Note that when two tracklets are Discordant (have 
detections in common) if one of the two is included in one confirmed identifi- 
cation then the other can be considered used. In this way the set of tracklets 
is reduced and many false tracklets are discarded. In the opposition simulation 
there are 168, 122 leftover tracklets of which 5, 363 are false: a reduction by 
74.3% of the original dataset and a reduction by 79.4% in the number of false 
tracklets. At the sweet spots the corresponding numbers are 232, 101 leftover 
tracklets of which 17, 033 are false: a reduction by 66.6% and 71.2% respectively. 

The leftover tracklets can be used as input to another iteration which could 
use the same algorithms as the first one (maybe with different controls and 
options) or could use very different methods; because of the reduced number 

^^This incompleteness was in fact a problem of interface between two processing steps: 
complete identifications could have been found by the track composition algorithms. 
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density of tracklets, Accuracy should be less of a problem. One possibility is 
to iterate the same procedure: generate tracks starting from the leftover track- 
lets with the same algorithm |Kubica et al. 2007] but using different options; 
then run again the orbit determination, possibly with different options, and the 
normalization. Another possibility is to use for a second iteration a completely 
independent algorithm to find identifications. Such an algorithm has already 
been proposed and tested on large simulations |Milani et al. 2005b| and on small 
real data sets jBoattini et al. 2007] . 

A full discussion of the iteration strategies to be used for the next gener- 
ation surveys is beyond the scope of this paper. However, to show that the 
normalized Efficiency values shown in Table [7] are not to be considered a prob- 
lem, we have run an improved version of the recursive attribution algorithms of 
[Milani et al. 2005bj on the leftover trackletQ. The results for the opposition 
simulation are as follows: we have been able to recover 75.4% of the lost objects 
(85.3% of the lost NEOs), thus bringing the overall Efficiency to 99.0% (97.1% 
for NEOs); for the sweet spots the values are 75.0% recovery (85.2% for NEOs) 
and 99.4% overall Efficiency (97.1% for NEOs). The false identifications remain 
very few even including the second iteration results: only 86 (0.06%) at oppo- 
sition and only 29 (0.02%) at the sweet spots. The same procedure allows also 
to compute normalized identifications for 2-nighters, with Efficiency 83.4% and 
only 2.1% false identifications at opposition. The corresponding values in the 
sweet spots are 89.2% Efficiency and 1.3% false identifications. 

The best way to assess Goodness of the results is to try to perform a sim- 
ulation mimicking as much as possible the way in which the results from one 
lunation would be used in the next one. Thus we have made two complete sim- 
ulations at opposition for two consecutive lunations. Having obtained 3-nighter 
identifications for the first month we have attempted to attribute to each one 
of them some tracklets in the other month. Using the same quality controls 
as for orbit determination within the same lunation this procedure was 99.6% 
efficient for objects with 1 tracklct in the second lunation, 99.7% efficient for 
objects with 2 tracklets and 99.9% efficient for objects with 3 tracklets. There 
were no NEOs among the few cases of either failed or incomplete attribution. 

7 Conclusions and Future Work 

The purpose of this paper was to identify efficient algorithms to compute pre- 
liminary and least squares orbits given a track (or proposed identification). 

We have found suitable algorithms by revising the classical preliminary orbit 
methods. The most important improvements are provisions to keep alternate 
solutions under control. The existence of double solutions was known since long 
time and we have shown that even triple solutions can occur. Still there is no 
reason this should impair the performance of the orbit determination algorithms, 
provided these cases are handled with due care. 

^^To better control the false identifications wc have used even tighter quality controls. 
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For the differential corrections stage, leading from preliminary orbits to least 
squares ones, we have adopted algorithms innovative with respect to the clas- 
sical ones but already available from our previous work. If properly combined 
with a good control logic they significantly improve the efficiency of differen- 
tial corrections even when the preliminary orbits are not close to the nominal 
solution (e.g., because of low curvature). 

The third stage of orbit determination is the quality control of the results 
imposed by applying statistical criteria to the residuals of the least squares fit. 
When there is just one object (or only a small number) under consideration at 
a time this stage may appear unnecessary. However, with the high sky-plane 
detection density to be expected with the next generation surveys this stage 
turns out to be the most critical one. When the number density of observations 
(per unit area on the sky) is large tracklets belonging to different objects may 
be incorrectly identified. To clean up the output from these false identifications 
is not easy. We have found the method of normalization to be very effective 
in removing false identifications, but unavoidably some true identifications are 
sacrificed to remove the discordant false ones. The critical point is to select 
options and details of the algorithms in such a way that the number of false 
identifications is kept to a very low value but few true ones are lost. 

Although our mathematically rigorous theoretical results do not need con- 
firmation it has been useful to test their practical performance on simulations 
of the next generation surveys. In this way we have shown that orbits can be 
computed even for the most difficult classes of orbits. We have also shown, 
with full density simulations including an overwhelming majority of MBA, that 
the large number of objects observed does not result in a "false identification 
catastrophe" . On the contrary, a large number density is compatible with a 
low number of lost objects provided the quality control on the residuals is tight 
enough and the sequence of algorithms is suitably chosen. 

The performance of a procedure for identification and orbit determination 
critically depends upon the performance of the individual algorithms and upon 
the pipeline design - the sequence of algorithms operating one upon the output 
of another. We have used the algorithms from [Kubica et al. 2007] as the first 
step, followed by the algorithms of this paper as the second step. We have 
mentioned the possibility of using the algorithms of [Milani et al. 2005b| as the 
third step. Even more complicated pipelines can be conceived and may be 
superior in performance. However, a detailed discussion of pipeline design is 
beyond the scope of this paper and will be the subject of future work. 

Another subject of future work is the definition of the procedure to combine 
the results from different lunations of a large survey. We have done a small test 
with a second lunation by using the attribution algorithm of [Milani et al. 2001] . 
Under other conditions, when a survey has been operating for more than one 
year, other algorithms such as the ones of [Milani et al. 200"0] and a variant of 
the methods of [Kubica et al. 2007] may become necessary. 

In the three papers [Milani et al. 2005b] , [Kubica et al. 2007] and the present 
one we have defined a number of algorithms to be used to process astrometric 
data of Solar System objects when the number density will be much larger than 
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that currently observed. This wiU very soon be the case with the next genera- 
tion surveys including Pan-STARRS and LSST. Such algorithm definition work 
is a necessary step to exploit their superior survey performance and provide 
identifications and orbits for most observed objects. In the near future we will 
need to handle real data which, of course, will contain unpredicted problems 
and present a new, formidable challenge. 
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